function idat = fsb_normalize_baseline(idat,sandbox)

% FSB :  Trial Baseline Normalization
%
% EXAMPLE:
% idat = fsb_normalize_baseline(idat,sandbox)
%
% INPUT:
% idat:         4-D image data
% sandbox:      sandbox experiment struct
%
% OUTPUT:
% idat:         Baseline corrected 4D image data
% 
% CALLED by:
% FSB.m
%
% NOTES:
% Due to datatype issues, using this function might result in artificial
% activations. Check for activations outside of the brain
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License: GNU GPL, no express or implied warranties
% 
%$ Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Normalizing trials with own baseline intensity...');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Transform int16 to single for further processing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
idat = single(idat);
idat(idat==0)=NaN;
idat(idat<=0)=NaN;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Define baseline period
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
blength=zeros(1,max(sandbox.intrial(:,3)));
vxcorr = zeros(1,max(sandbox.intrial(:,3)));
scanind = cell(1,max(sandbox.intrial(:,3)));
bsl_ind = scanind;
for x = 1:max(sandbox.intrial(:,3));

    scanind{x} = find (sandbox.intrial(:,3)==x);
    bsl_trial = sum(sandbox.hemodynamics(scanind{x},1:4),2); % Doing alternative baseline calculation, stimulation regressors only
    bsl_max = max(bsl_trial);
    ind_t = find(bsl_trial<(bsl_max/20));
    if ~isempty(ind_t);
        bsl_ind{x} = scanind{x}(ind_t);
    end

end

bsl_index = vertcat(bsl_ind{:});
bsl_idat = idat(:,:,:,bsl_index);
bsl_mean = nanmean(bsl_idat,4);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Do actual normalization
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

for x = 1:size(scanind,2);

    waitbar(x/max(sandbox.intrial(:,3)));
    idat_red = idat(:,:,:,scanind{x}); % Gives volumes in trials

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Looks for baseline volumes by looking for values that are below
    % 1/20th of maximum response
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if ~isempty(bsl_ind{x});
        bsl_ind2 = mean(idat(:,:,:,bsl_ind{x}'),4);
        voxel_corr = bsl_mean./bsl_ind2; % correction factor voxel by voxel
        voxel_corr(voxel_corr==inf)=1;

        voxel_corr_4D = repmat(voxel_corr,[1 1 1 size(idat_red,4)]);
        idat_red = idat_red.* voxel_corr_4D;
        idat(:,:,:,scanind{x}) = idat_red;
        vxcorr(x) = nanmean(voxel_corr(:));
        blength(x)= size(bsl_ind,1);
    end
end

vxcorr = vxcorr/nanmean(idat(:))*100;
disp(['Mean baseline length: ' num2str(mean(blength))]);
disp(['mean voxel corr: ' num2str(nanmean(vxcorr))]);
close(h);
figure; bar(vxcorr);
set(gcf,'Name','Trial baseline intensity correction');
xlabel('Trial');
ylabel('Average signal deviation %');
disp(' ')
end
